Ground state description of a single vortex in an atomic Fermi gas: From BCS to Bose-Einstein 

condensation 
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We use a Bogoliubov-de Gennes (BdG) formulation to describe a single vortex in a neutral fermionic gas. It is 
presumed that the attractive pairing interaction can be arbitrarily tuned to exhibit a crossover from BCS to Bose- 
Einstein condensation. Our starting point is the BCS-Leggett mean field ground state for which a BdG approach 
is microscopically justified. At strong coupling, we demonstrate that this approach is analytically equivalent to 
the Gross-Pitaevskii description of vortices in true bosonic systems. We analyze the sizable density depletion 
found for the unitary regime and relate it to the presence of unoccupied (positive energy) quasi-bound states at 
the core center 



PACS numbers: 03.75.Hh, 03.75.Ss, 74.20.-z 

One of the most exciting developments in atomic and con- 
densed matter physics has been the observation of superfluid- 
ity in trapped fermionic systems | llQ-IHBl- In these systems, 
the presence of a Feshbach resonance provides a means of 
tuning the attractive pairing interaction with applied magnetic 
field. In this way the system undergoes a continuous evolution 
from BCS to Bose-Einstein condensed (BEC) superfluidity. 

The most conclusive demonstration of the superfluid phase 
has been the experimental observation of vortices |5]. Partic- 
ularly interesting from a theoretical viewpoint is the way vor- 
tices evolve from BCS to BEC. This evolution is associated, 
not just with a decrease in vortex size but with a complete rear- 
rangement of the fermionic states which make up the core. As 
a result, there is a continuous evolution of the particle density 
within a vortex, thereby affecting the visibility of vortices in 
the laboratory. In this paper we discuss the behavior of a (sin- 
gle) vortex as the system crosses from BCS to BEC. Our work 
is based on simplest BCS-like ground state first introduced 
by Leggett |6] and Eagles |7] to treat BCS-BEC crossover. 
With this choice of ground state inhomogeneity effects are 
readily incorporated as in generalized Bogoliubov-de Gennes 
(BdG) theory. Here we demonstrate analytically that the BdG 
strong coupling description of the T — Q vortex state coin- 
cides with the usual Gross-Pitaevskii (GP) treatment of vor- 
tices in bosonic superfluids. A fermionic theory based on BdG 
is, thus, very inclusive, and within this approach one expects 
a smooth evolution of vortices from the BCS to BEC limit as 
the statistics effectively change from fermionic to bosonic. 

Previous studies of vortices in these fermionic superfluids 
addressed the BCS limit at T = |8] and T w |9]. There 
is also work 1 10] on the T = strict unitary case where a 
BdG approach was used with Hartree-Fock contributions in- 
cluded. In the present work, by contrast, we discuss the entire 
crossover regime and, importantly, present a detailed analysis 
of the energy and spatial structure within the core and how it 
evolves from BCS to BEC. A very different path integral ap- 
proach was introduced in Ref. Illll to address vortices with 
BCS-BEC crossover, but here the authors note that density 
depletion effects appear to be unphysically large in the BCS 
regime. Our analytical approach builds heavily on previous 
work 1 12] which showed a general connection between GP 
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theory and BdG. From this one can conclude that a general- 
ized BCS theory ]6] treats the bosonic degrees of freedom at 
the same level as GP theory. Different ground states can be 
contemplated, (with incomplete condensation, say) but they 
will not be compatible with BdG theory. In a similar way, 
once T 7^ one has to incorporate noncondensed pairs, and 
associated pseudogap physics USUI which are not present in a 
finite temperature BdG theory. 

For the mostpart, BdG approaches require detailed numer- 
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ical solution ]8, 14, 15, 16], so it is particularly useful to 
have analytical tools in the BEC limit. We present this non- 
numerical description first. Our general self consistent equa- 
tions 11711 are 
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-^V^+Ve^tir), A(r) is the T = gap function 
which is importantly the same as the T = order parameter, 
Asc(r), Vea;t(r) is the external potential associated with the 
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trap, and we choose h = 1 with J dr{ 



for all energy levels n. The difference between the present 
approach and the usual BdG applications of superconductivity 
is that here the fermionic chemical potential /i must be self 
consistently determined, as the attractive coupling constant is 
varied. 



We use a Green's function formulation 11411 to write the 
zero-temperature free energy Eq = {H — nN), (where N 
is the number operator) in the form 
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Here (i=l,2,3) are the Pauli matrices. The elements of G 
corresponding to the normal and anomalous channels can be 
further expressed as coupled integral equations in terms of the 
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noninteracting Green's function G\i 

Gn(r,r';c.) = G^"^ (r, r'; a;) + J dr2[G^°,\r,r2;cj) 

X A(r2)G2i(r2,r';^)], (3) 

G2i(r, r'; ^) = " / '^r2Gl°^ (r2, r; -w) A*(r2)Gii(r2, r'; 

(4) 

where the gap and number equation are given by 
A*(r) = -Vj^^§^G2i{r,r;cu) and n(r) = 
2 f^e*"''Gii(r,r;w). The coupHng constant ~V < 
for the attractive inter-fermion contact interaction is 
parametrized in terms of the s-wave scattering length a p with 
m/A-KaF ^ -l/V + J2k m/P. 

In the strong pairing limit a^? is small and positive and 
A/|/i| is small. To derive the ground state energy Eq from 
Eq.(|3, we expand the Green's function G in the gap equation 
in terms llT2ll of G^ii . Including terms up to fourth order in 
A|, it follows that 

Eo[A] = I dr{^i^-ao(r)|A(r)p 

+ i6o(r)|VA(r)p_ico(r)|A(r)|4} (5) 

where ao(r) ^ + ^[Mi3 " 2V^e.t(r)], 6o(r) =^ 
maF/ldiT, and co(r) ~ —m^ap/Wn. Here = 2/z + eg 
is the effective "bare" bosonic chemical potential, and £o = 
{ma^)^^ is the binding energy of the composite boson, with 
I /is I ^ £o- It is assumed that Vext is slowly varying 1121 so 
thatG(°)(r,r';c^,) ~ G'([\r-r'; fi-[Vext{r) + Vext{r')]/2). 
Similarly, A(r) is assumed to vary slowly on the scale of ap. 
As a consequence of these assumptions, and for the purposes 
of this paper, (which ultimately focuses on a single vortex), 
trap effects are not particularly relevant. 

It should be stressed that this expansion is similar to 
Gor'kov's derivation of Ginzburg-Landau (GL) theory in the 
BCS limit at T « Tc, albeit here we consider strong coupling 
and T = 0. As in conventional superconductors |18], mini- 
mizing the energy £'o[^'] with respect to ^E", one obtains 

Eo[-^] = f dr 15-14 + ^ f dl ^*n.V*|s. (6) 

niB J 2m_B J 

where we have identified the condensate wave function ^I* (r) 
as \Jrf?ap l^-K A(r) and h is the unit vector normal to the 
surface. The second term is a surface term, which vanishes 
for an infinite system. In a neutral superfluid, however, the 
energy has to be calculated within a finite region of radius R 
around the vortex core, to avoid divergences, so this surface 
term cannot be neglected. 

Equivalently, the zero-temperature energy can be written in 
a more conventional form as 

E^{^\ = y'dr{^|V^'(r)p + 2Fe.t(r)|vi/(r)p 

+ i[/o|vI/(r)r~^5|*(r)p}, (7) 



where [/q = AiraB /mB, tub = 2m and = 2aF are the 
mass and the scattering length of the composite boson. 

Importantly, this expression has the same form as the T — Q 
energy of a gas of weakly-interacting bosons {l^ associated 
with GP theory. It should be noted that here this GP theory is 
written in terms of the grand canonical representation where 
the bosonic chemical potential (rather than the number of par- 
ticles N) is held fixed. Minimizing the zero-temperature en- 
ergy Eq (via 5Eq/5'^* = 0), leads to the well known GP 
equation: 

-77^V2*(r) + 2Vext{v) + (7o|*(r)|2*(r) = mb^-W. 
2mB 

(8) 

We emphasize that this BdG analysis has, in effect, derived 
GP theory from a fermionic starting point. The presence of a 
Hartree term in the BdG equations will destroy the simple an- 
alytic arguments presented here. Nevertheless, a Hartree con- 
tribution for the composite bosons is found here of the form 
J/ol^'l^^ — ^ UqUb'^', (where ub is the density of bosons). 
This is, of course, unrelated to the Hartree term of the origi- 
nal fermions, which is absent in Eq.Q, as is consistent with 
the usual BCS-Leggett mean field theory |6]. The inclusion 
of Hartree terms 1 8 ] in the vortex problem has been accounted 
for in the literature, at weak |8| and strict unitary coupling 
1 10], but they do not appear to lead to important differences. 

Expanding the zero-temperature energy Eq to next order in 
A, a term /drdo(i")|A(r)|'^ with dp (r) - -5m''^a^/2567r 
appears in. This term contributes a term (73|5'(r)|4vl/(r) with 

— — 157r^a^/4mB in Eq. it introduces the appropriate 
analogue in Eq. Q as well. This three-body correction to the 
usual GP equation (g^) Jl^l represents an effective attraction. 
In the composite-boson system, it provides a first-step correc- 
tion of the BEC limit en route to the fermionic or BCS end 
point. 

In the presence of a single vortex, the wave function 
can be written as ^'(r) = /(r)e^*^. We introduce the 
BEC correlation length ^bec in the strong-pairing limit as 
{2mB^BEcy^ = A*s = UafS, with /o ^ f{r oo). We 
rescale r — ^bec • x and fjr) — fo ■ X(x) and apply stan- 
dard boundary conditions I191 120I1 . The results for \['(r) in 
these units are plotted as the solid lines in Fig. [2 As shown 
in Fig.[^, (and consistent with earlier results in the literature 
20]), the wavefunction rises smoothly from zero at the 
center of the core to its full magnitude at infinity on a length 
scale of ^BEC- 

An important feature of this figure should be noted. In this 
BEC limit of the BdG equations, the wavefunction is smooth. 
This behavior, which is in contrast to the BCS limit, is a con- 
sequence of the absence of localized fermionic bound states in 
the core region. In the weak coupling limit, as first noted by 
Caroli et al i2ll . these bound states are associated with energy 
eigenvalues En < Aqo- Here Aqo is the value of gap function 
in the bulk, away from the core. The gap A(r) provides an ef- 
fective potential well for the quasiparticles around the vortex 
core. From Eq. Q we have En > — /i. Since — /i ^ Aoo, 
these bound states are necessarily absent in the strong pairing 
limit . 

The first appearance of fermionic properties in our 
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Figure 1: (a) Numerical solution of GP equation in the BEC limit with (g'^ = 0.1, dashed lines) and without (solid lines) the three-body 
ga-term, (b) the corresponding normalized particle density n{x)/noo as a function of x, and (c) the zero-temperature vortex energy cost Ev as 
a function of R/^bec- In (c), the difference is shown as the red dot-dashed curve . Here Uoo = n(oo). 




Figure 2: Local fermionic density of states N{E,r) as function of E for BCS {kpa^ — 1), unitary and BEC {kpaKiV) cases, at the center 
r = Q (black dashed curves) and radius r = 25/fcF (red solid curves) of the vortex core. The bulk value of the gap Aoo is 0.21, 0.68, I.TiEf, 
respectively. In the BEC case, fj, ~ —Q.8Ef. Here Ef is the noninteracting Fermi energy at the trap center. 



composite-boson system is via the ga-term in the GP equa- 
tion and the corresponding term in the zero-temperature en- 
ergy E^. The effects of this addition are plotted as dashed 
lines in Fig.[2for (rescaled) g'^ = ~93fo /Uq = 0.1. This 93 
contribution represents an attractive interaction, and, as shown 
in the figure, leads to a slight increase in the core size. 

One can similarly compute the particle density ns (r) asso- 
ciated with composite bosons which is simply related to the 
wave function as nsir) = |5'(r)p = n(r)/2, where n is the 
density of fermions. The density n{x) is plotted in Fig.^b), 
normalized at the bulk value = n{oo), as a function of 
X. As expected the particle density is strictly zero at the core 
center in the BEC limit, where there is complete depletion. 

The energy cost of a single vortex can also be calcu- 
lated from Eq. and the result is plotted in Fig. 1(c). 
The energy cost per unit length is given by E^, = 



xdx, where R is 



ms Jo [\dx) ~'"^~'^2V^ 

a cutoff needed to regularize a calculation of the vortex core 
energy. In Fig. ^c) the solid line indicates E^ as a func- 
tion of R/^BEC- The shape of the curve at the region 
R/£,BEC > 2 can be fitted to the usual functional form 
Ey oc IniDcR/^BEc), where Dq = 1.48. The dashed line 
in Fig.^c) presents results for Ey in the presence of the three 
body term, where we take g'^ — 0.1. This correction (red 
dot-dashed curve) lowers the vortex energy, as shown in the 
figure, and it approaches an asymptote as R/^bec ~> 00. 



To understand the details of the core structure we turn now 
to numerical solutions of the BdG equation. We build on our 
analytical analysis at strong coupling to provide a check for 
our numerical algorithms. Here the physical coupling strength 
is controlled by the parameter 1/kpa. We have verified that 
changes in our high cutoff energy and associated coupling 
constant V do not affect the numerical results. Our numerical 
method is very similar to that in Ref. |22J. The chemical po- 
tential /i is approximated by the homogeneous solution, since 
the vortex core only occupies a small portion of the entire sys- 
tem. Here we begin with a study of the localized (fermionic) 
density of states (LDOS), N{E,r), within the core region. 
There is considerable interest in the literature in the behavior 
of the LDOS for high Tc |23] as well as low Tc supercon- 
ductors 1 15], since this quantity is accessible through scan- 
ning tunneling microscopy measurements. N{E, r) is given 
by EnKS{E - En) + vlS{E + En)]. We ignore, for nu- 
merical simplicity, dependencies of the wave functions on the 
cylindrical variable z, since these do not lead to qualitative 
effects. Integrating N{E, r) over _E < reflects the particle 
density distribution n{T) inside the core. Thus, this quantity 
provides a means of understanding the density depletion, or 
lack thereof, inside the core. It is essential for arriving at a 
deeper understanding of the core region and structure. 

In Fig.|2]we plot N{E, r) inside the core, as a function of 
energy E for r — and r — 25/kp. The three panels corre- 
spond to BCS (the noninteracting wavevector at the trap center 
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Figure 3: Gap function A (black dashied curves) and particle density 
n (red solid curves) as a function of r for BCS {kpa — —1) and 
unitary cases. All quantities are normalized by their bulk values. 
The inserts in the left panels shows detailed behavior of the gap near 
the center of the vortex core. 

kpa = —1), unitary and BEC (kpa = 1). Rather than show- 
ing resuhs for the strict BCS and BEC regimes, these plots 
represent a physically accessible range of magnetic fields. Be- 
cause N{E, r) is a fundamentally fermionic quantity, this in- 
formation is lost from the analytical analysis which leads to 
a Gross-Pitaevskii transcription of the BEC limit. In the BCS 
case, the r — peak atE — reflects a bound fermionic state. 
This state together with the continuum of scattering states is 
responsible for the fact that there is no core density depletion. 
For the unitary regime, the lower energy peak in Fig.|2]arises 
from scattering states and appears at energies near the bulk 
gap Aoo- This quasi-bound state has energy close to the scat- 
tering state continuum and is reflected in the LDOS by a slight 
movement of the BCS central peak to the right. The energy 
integral of this feature is important in determining the finite 
particle density n(r) at the core center. The peak at positive 



energy is a reflection of a quasi-bound state. This unoccupied 
quasi-bound state effectively depletes the spectral weight for 
E < 0, and therefore leads to the density depletion within 
the core. By the BEC limit all remnants of fermionic states 
have disappeared, until one provides energy large enough to 
break the pairs. It can be seen from the figure that in all three 
cases at sufficiently large distances from the core center the 
fermionic density of states assumes the bulk value. 

In Fig. 3 we plot the position dependent order parameter 
A(r) along with the particle density distribution n(r) for the 
unitary and BCS (kpa ~ —1) cases shown in the previous 
figure. This BCS-like case still has a reasonably large bulk 
gap, so there a non-negligible depletion at the core. This is 
in contrast to arbitrarily weak coupling, where the depletion 
vanishes. The small r oscillations shown here in both A(r) 
and n(r) reflect the presence of a true bound state, as in earlier 
work ■ The oscillations at large r are an effect of the finite 
system size. 

Importantly, in the unitary case, the particle density at the 
core center is substantially lower than in the BCS case. This is 
a consequence of the reduced spectral weight seen in the lower 
peak in the middle panel of Fig. 2. Our results are within 20% 
of those obtained in Ref. fioll . In this earlier work a Hartree- 
Fock correction was applied to the BdG equations which was 
argued |24] to be the source of the density depletion. Here, 
we interpret this depletion differently as associated with the 
behavior of the core excitation spectra in conjunction with the 
reduced chemical potential (/i < Ep)- 

We thank N. Nygaard for sharing his thesis. This work was 
supported by NSF-MRSEC Grant No. DMR-0213745. 

Note After this work was complete we learned of a related 
calculation by Machida and Koyama (Phys. Rev. Lett. 94, 
140401 (2005)) which attributed the density depletion within 
the vortex core at unitarity to closed-channel bosons. 
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